/* visFlowTopology.cpp
*/
#include "visFlowTopology.h"

namespace visualization
{
void visFlowTopology::determineCriticalPoint()
{
//	cout << "Input the interpolation neighbour number, it can be 4 or 5..." << endl;
//	cin >> wunai;

	int k = 0;
	double pos[3];
	// Scan the point set.
	for(size_t i = 1; i <= dataField->getPntNumberOfModel(); i++)
	{
		if(isZeroPoint(i, pos))
		{
			k++;
			no_zero_point = 0;
			if(!computeJacobian(pos))	{ continue;	}
			if(matrixIsZero(jacobian)) { continue; }
			vis2by2RealMatrixEE(jacobian, eigenvector, eigenvalue);
//			if(!calculator(jacobian, eigenvalue, eigenvector)) { exit(-1); }
			if(classify(pos))
			{
				/**Since only the newly-added-saddle( in function classify() )'s infomation is needed,
				 * so the accessed item is "saddles.size() - 1".
				 */
				for(int j = 0; j < dim * dim; j++)
				{
					//In fact, I'm not sure what if the imaginary part of the eigenvector is not zero.
					saddles[saddles.size() - 1].eigenvec[j] = eigenvector[j].real();
				}
				// Normalize it.
				normalizeVec( saddles[saddles.size() - 1].eigenvec );
			}
		}
	}
}

void visFlowTopology::preProcessSaddles()
{
	groupSaddles();
	removePesudo();
	adjustSaddlePosition();
}

void visFlowTopology::groupSaddles(void)
{
	// Within this distance, the two saddles are consider to potentially as the same one saddle.
	const double distance = 5 * radius;
	const double z = 0.6;
	for(unsigned int i = 0; i < saddles.size(); i++)
	{
		// Indexes of saddles in variable saddles that are potentially the same one.
		vector<int> oneGroup;
		for(unsigned int j = 0; j < saddles.size(); j++)
		{
			/* 1. compute distance differences.
			 * 2. compare eigenvector differences.
			 *    Note that here the eigenvectors have been normalized when the program is running.
			 *    if they're the same one saddle, push the index j to the back of oneGroup.
			 */

			//1.
			double dis = sqrt( pow( saddles[j].x - saddles[i].x, 2 ) + 
									  		 pow( saddles[j].y - saddles[i].y, 2 )
											 );

			//2.
			//2.1 compute according to the original order of the two eigenvactors.
			/*    /        \
			 *    |  0   2 |
			 *    |  1   3 |
			 *    \        /
			 */
			double dotProduct[2];
			dotProduct[0]  = saddles[i].eigenvec[0] * saddles[j].eigenvec[0];
			dotProduct[0] += saddles[i].eigenvec[1] * saddles[j].eigenvec[1];
			dotProduct[1]  = saddles[i].eigenvec[2] * saddles[j].eigenvec[2];
			dotProduct[1] += saddles[i].eigenvec[3] * saddles[j].eigenvec[3];
			if( abs(dotProduct[0]) > z && abs(dotProduct[1]) > z && dis < distance )
			{
				if(dotProduct[0] < 0)
				{
					saddles[j].eigenvec[0] *= -1;
					saddles[j].eigenvec[1] *= -1;
				}
				if(dotProduct[1] < 0)
				{
					saddles[j].eigenvec[2] *= -1;
					saddles[j].eigenvec[3] *= -1;
				}
				oneGroup.push_back(j);
				continue;
			}
			        
			//2.2 compute according to the opposite order of the two eigenvectors.
			dotProduct[0]  = saddles[i].eigenvec[0] * saddles[j].eigenvec[2];
			dotProduct[0] += saddles[i].eigenvec[1] * saddles[j].eigenvec[3];
			dotProduct[1]  = saddles[i].eigenvec[2] * saddles[j].eigenvec[0];
			dotProduct[1] += saddles[i].eigenvec[3] * saddles[j].eigenvec[1];

			if( abs(dotProduct[0]) > z && abs(dotProduct[1]) > z && dis < distance )
			{
				if(dotProduct[0] < 0)
				{
					saddles[j].eigenvec[2] *= -1;
					saddles[j].eigenvec[3] *= -1;
				}
				if(dotProduct[1] < 0)
				{
					saddles[j].eigenvec[0] *= -1;
					saddles[j].eigenvec[1] *= -1;
				}
				// Exchange the eigenvectors of jth saddle.
				double tmp[2] = { saddles[j].eigenvec[0], saddles[j].eigenvec[1] };
				saddles[j].eigenvec[0] = saddles[j].eigenvec[2];
				saddles[j].eigenvec[1] = saddles[j].eigenvec[3];

				saddles[j].eigenvec[2] = tmp[0];
				saddles[j].eigenvec[3] = tmp[1];
				// Push the j into oneGroup.
				oneGroup.push_back(j);
				continue;
			}
		}//for
		if(1 < oneGroup.size())
		{
			unsigned int size = oneGroup.size();
			// Save the average location and eigenvectors in the ith of the variable saddles.
			// 1. save the total to the ith of it;
			// 2. compute the average.
			// 3. erase unusefull items.

			// 1
			for(unsigned int k = 1; k < size; k++)
			{
				saddles[i].x += saddles[ oneGroup[k] ].x;
				saddles[i].y += saddles[ oneGroup[k] ].y;
				saddles[i].z += saddles[ oneGroup[k] ].z;

				for(int d = 0; d < dim * dim; d++)
				{
					saddles[i].eigenvec[d] += saddles[ oneGroup[k] ].eigenvec[d];
				}
			}
			// 2
			saddles[i].x /= size;
			saddles[i].y /= size;
			saddles[i].z /= size;
			for(int d = 0; d < dim * dim; d++) {saddles[i].eigenvec[d] /= size;}
			// 3
			for(unsigned int k = 1; k < size; k++)
			{
				int index = oneGroup[k] - (k - 1);
				saddles.erase( saddles.begin() + index );
			}
		}//for
	}//for
}
void visFlowTopology::removePesudo(void)
{
	// Remove pesudo-saddle
	/* Four integration starting points.
	 *             d2F
	 *             /|\
	 *              |
	 *    d1B --->saddle<--- d1F --'pnt'---
	 *              |
	 *             \|/
	 *              d2B
	 */

	cout << "Input threshold for topology(projection)... " << endl;
//const double diff = 0.75; // 0.75for data 2
//const double diff = 0.5; // 0.5for data 19
	double diff;
	cin >> diff;
	for(size_t i = 0; i < saddles.size(); i++)
	{
		double pos[3] = {saddles[i].x, saddles[i].y, saddles[i].z};
		double prj[4];
		projection(pos, i, 1, prj);
				if( abs(prj[0]) < diff || abs(prj[2]) < diff || prj[0] * prj[1] > 0 
		   ||
				abs(prj[1]) < diff || abs(prj[3]) < diff || prj[2] * prj[3] > 0
		  )
		{
			saddles.erase(saddles.begin() + i);
			i--;
		}
	}
}

void visFlowTopology::adjustSaddlePosition(void)
{
	const double exactPos = 0.95;
	const double shiftStepLen = 0.01;
	for(size_t i = 0; i < saddles.size(); i++)
	{
		// 1. Imagic that shift the saddle along each of its eigenvectors by a distance;
		// 2. Compute the smmallest projection;
		// 3. Really shift the saddles to the imagic position;
		// 4. Repeat 1~3 untill the projection is less than exactPos(the threshold).
		double originalPos[3] = { saddles[i].x, saddles[i].y, saddles[i].z };
		double maxPrj = DBL_MAX;
		do
		{
			double shift1[3] = { originalPos[0]+shiftStepLen, originalPos[1], 0};
			double shift2[3] = { originalPos[0]-shiftStepLen, originalPos[1], 0};
			double shift3[3] = { originalPos[0], originalPos[1]+shiftStepLen, 0};
			double shift4[3] = { originalPos[0], originalPos[1]-shiftStepLen, 0};
//			double shift1[3] = { originalPos[0]+saddles[i].eigenvec[0] * shiftStepLen,
//													 originalPos[1]+saddles[i].eigenvec[1] * shiftStepLen, 
//													 0
//												 };
//			double shift2[3] = { originalPos[0]-saddles[i].eigenvec[0] * shiftStepLen,
//													 originalPos[1]-saddles[i].eigenvec[1] * shiftStepLen, 
//												   0
//												 };
//			double shift3[3] = { originalPos[0]+saddles[i].eigenvec[2] * shiftStepLen,
//												   originalPos[1]+saddles[i].eigenvec[3] * shiftStepLen, 
//													 0
//												  };
//			double shift4[3] = { originalPos[0]-saddles[i].eigenvec[2] * shiftStepLen,
//			                     originalPos[1]-saddles[i].eigenvec[3] * shiftStepLen, 
//													 0
//												 };
//
			// projections of four directions.
			double origPrj[4], prj1[4], prj2[4], prj3[4], prj4[4];
			projection(originalPos, i, diavation, origPrj);
			projection(shift1, i, diavation, prj1);
			projection(shift2, i, diavation, prj2);
			projection(shift3, i, diavation, prj3);
			projection(shift4, i, diavation, prj4);
			double prjSum[5] = {0, 0, 0, 0, 0};
			
			for(size_t i = 0; i < 4; i++)
			{
				prjSum[0] += abs(origPrj[i]); 
				prjSum[1] += abs(prj1[i]);
				prjSum[2] += abs(prj2[i]);
				prjSum[3] += abs(prj3[i]);
				prjSum[4] += abs(prj4[i]);
			}
			size_t shiftTo = 0;
			for(size_t i = 1; i < 5; i++)
			{
				shiftTo = prjSum[shiftTo] > prjSum[i] ? shiftTo : i;
			}
			
			switch(shiftTo)
			{
				case 0:
					maxPrj = 1;
					break;
				case 1:
					maxPrj = prjSum[1]/4;
					for(size_t j=0;j<3;j++)originalPos[j]=shift1[j];
					break;
				case 2:
					maxPrj = prjSum[2]/4;
					for(size_t j=0;j<3;j++)originalPos[j]=shift2[j];
					break;
				case 3:
					maxPrj = prjSum[3]/4;
					for(size_t j=0;j<3;j++)originalPos[j]=shift3[j];
					break;
				case 4:
					maxPrj = prjSum[4]/4;
					for(size_t j=0;j<3;j++)originalPos[j]=shift4[j];
					break;
			}
		}while(maxPrj < exactPos);
		cout << maxPrj << endl;
		saddles[i].x = originalPos[0];
		saddles[i].y = originalPos[1];
		saddles[i].z = originalPos[2];
	}
}

void visFlowTopology::projection(const double pos[3], const size_t sdl, const double dist, double prj[4])
{
	double xDiv = saddles[sdl].eigenvec[0] * dist;
	double yDiv = saddles[sdl].eigenvec[1] * dist;
	double d1F[2] = { pos[0] + xDiv, pos[1] + yDiv };
	double d1B[2] = { pos[0] - xDiv, pos[1] - yDiv };
	xDiv = saddles[sdl].eigenvec[2] * dist;
	yDiv = saddles[sdl].eigenvec[3] * dist;
	double d2F[2] = { pos[0] + xDiv, pos[1] + yDiv };
	double d2B[2] = { pos[0] - xDiv, pos[1] - yDiv };

	double pnt1[3] = {d1F[0], d1F[1], 0};
	vector<double> *valueOnpnt1 = new vector<double>;
	if( !dataField->interpolatedValue(pnt1, 3, valueOnpnt1, 5, 1) )
	{
		cout << "Sorry! It should fail to interpolate... at line " 
				 << __LINE__ << "	of file " << __FILE__ << endl;
		exit(-1);
	}
	double pnt3[3] = {d1B[0], d1B[1], 0};
	vector<double> *valueOnpnt3 = new vector<double>;
	if( !dataField->interpolatedValue(pnt3, 3, valueOnpnt3, 5, 1) )
	{
		cout << "Sorry! It should fail to interpolate... at line " 
				 << __LINE__ << "	of file " << __FILE__ << endl;
		exit(-1);
	}
	double pnt2[3] = {d2F[0], d2F[1], 0};
	vector<double> *valueOnpnt2 = new vector<double>;
	if( !dataField->interpolatedValue(pnt2, 3, valueOnpnt2, 5, 1) )
	{
		cout << "Sorry! It should fail to interpolate... at line " 
				 << __LINE__ << "	of file " << __FILE__ << endl;
		exit(-1);
	}
	double pnt4[3] = {d2B[0], d2B[1], 0};
	vector<double> *valueOnpnt4 = new vector<double>;
	if( !dataField->interpolatedValue(pnt4, 3, valueOnpnt4, 5, 1) )
	{
		cout << "Sorry! It should fail to interpolate... at line " 
				 << __LINE__ << "	of file " << __FILE__ << endl;
		exit(-1);
	}
	// Contert it to a vector.
	pnt1[0] -= pos[0];
	pnt1[1] -= pos[1];
	pnt3[0] -= pos[0];
	pnt3[1] -= pos[1];

	pnt2[0] -= pos[0];
	pnt2[1] -= pos[1];
	pnt4[0] -= pos[0];
	pnt4[1] -= pos[1];

	pnt1[0] /= dist;
	pnt1[1] /= dist;
	pnt3[0] /= dist;
	pnt3[1] /= dist;

	pnt2[0] /= dist;
	pnt2[1] /= dist;
	pnt4[0] /= dist;
	pnt4[1] /= dist;

	double m1 = sqrt( pow((*valueOnpnt1)[0], 2) + pow((*valueOnpnt1)[1], 2) );
	double m3 = sqrt( pow((*valueOnpnt3)[0], 2) + pow((*valueOnpnt3)[1], 2) );
	double m2 = sqrt( pow((*valueOnpnt2)[0], 2) + pow((*valueOnpnt2)[1], 2) );
	double m4 = sqrt( pow((*valueOnpnt4)[0], 2) + pow((*valueOnpnt4)[1], 2) );

	valueOnpnt1[0][0] /= m1;
	valueOnpnt1[0][1] /= m1;
	valueOnpnt3[0][0] /= m3;
	valueOnpnt3[0][1] /= m3;

	valueOnpnt2[0][0] /= m2;
	valueOnpnt2[0][1] /= m2;
	valueOnpnt4[0][0] /= m4;
	valueOnpnt4[0][1] /= m4;

	// Calculate the dotproduct of vector 'pnt1' and 'valueOnpnt1'
	prj[0] = pnt1[0] * (*valueOnpnt1)[0] + pnt1[1] * (*valueOnpnt1)[1];
	prj[2] = pnt3[0] * (*valueOnpnt3)[0] + pnt3[1] * (*valueOnpnt3)[1];
	prj[1] = pnt2[0] * (*valueOnpnt2)[0] + pnt2[1] * (*valueOnpnt2)[1];
	prj[3] = pnt4[0] * (*valueOnpnt4)[0] + pnt4[1] * (*valueOnpnt4)[1];
}

void visFlowTopology::integrateCurves(intgResults *&forward, intgResults *&backward)
{
	size_t saddleNum = saddles.size();
	cout << "saddles            : " << saddleNum << endl;
	if(0 == saddleNum)
	{
		cout << endl << "There're no saddle points, stop!" << endl << endl;
		exit(-1);
	}

	// Here 4 means that for two dimensional data field, a saddle point has 4 related 
	// integration curves.
	ForwardSources = new Point[saddleNum * 2];
	BackwardSources = new Point[saddleNum * 2];
	for(size_t i = 0; i < saddleNum; i++)
	{
		/** Four integration starting points.
		 *            d2F
		 *            /|\
		 *             |
		 *   d1B --->saddle<--- d1F --'pnt'---
		 *             |
		 *            \|/
		 *             d2B
		 */
		double xDiv = saddles[i].eigenvec[0] * diavation;
		double yDiv = saddles[i].eigenvec[1] * diavation;
		double d1F[2] = { saddles[i].x + xDiv, saddles[i].y + yDiv };
		double d1B[2] = { saddles[i].x - xDiv, saddles[i].y - yDiv };
		xDiv = saddles[i].eigenvec[2] * diavation;
		yDiv = saddles[i].eigenvec[3] * diavation;
		double d2F[2] = { saddles[i].x + xDiv, saddles[i].y + yDiv };
		double d2B[2] = { saddles[i].x - xDiv, saddles[i].y - yDiv };
	
		double pnt[3] = {d1F[0], d1F[1], 0};
		vector<double> *valueOnpnt = new vector<double>[3];
		if( !dataField->interpolatedValue(pnt, 3, valueOnpnt, 5, 1) )
		{
			cout << "Sorry! It should fail to interpolate... at line " 
					 << __LINE__ << "	of file " << __FILE__ << endl;
			exit(-1);
		}
		// Contert it to a vector from saddle to original pnt.
		pnt[0] -= saddles[i].x;
		pnt[1] -= saddles[i].y;

		double m1 = sqrt( pow((*valueOnpnt)[0], 2) + pow((*valueOnpnt)[1], 2) );
		valueOnpnt[0][0] /= m1;
		valueOnpnt[0][1] /= m1;
		// Calculate the dotproduct of vector 'pnt' and 'valueOnpnt'
		double dp = pnt[0] * (*valueOnpnt)[0] + pnt[1] * (*valueOnpnt)[1];

		if(0 < dp)
		{
			ForwardSources[2*i].setCoordinates(d1F[0], d1F[1], 0);
			ForwardSources[2*i + 1].setCoordinates(d1B[0], d1B[1], 0);		
			BackwardSources[2*i].setCoordinates(d2F[0], d2F[1], 0);
			BackwardSources[2*i + 1].setCoordinates(d2B[0], d2B[1], 0);
		}
		else
		{
			ForwardSources[2*i].setCoordinates(d2F[0], d2F[1], 0);
			ForwardSources[2*i + 1].setCoordinates(d2B[0], d2B[1], 0);
			BackwardSources[2*i].setCoordinates(d1F[0], d1F[1], 0);
			BackwardSources[2*i + 1].setCoordinates(d1B[0], d1B[1], 0);
		}
		delete []valueOnpnt;
	}
	/* Before, the integration direction is just forward, so the result is somehow 
	 * odd. This is done because I believe that the direction should be one-trip.
	 * But, in fact, it should be done in the double directions.
	 */
	cout << endl << "It is now integrating..." << endl << endl;
	cout << "Input the parameter file for streamline generation ..." << endl;
	string pfile;
	cin >> pfile;

	integrationForwardCurves.parseParameter(pfile);
	integrationForwardCurves.setIntegrationForward();
	integrationForwardCurves.setSources(ForwardSources, saddleNum*2);
	integrationForwardCurves.setParameters(0.5, 2000);
  forward = integrationForwardCurves.integrate();
	
	integrationBackwardCurves.parseParameter(pfile);
	integrationBackwardCurves.setIntegrationBackward();
	integrationBackwardCurves.setSources(BackwardSources, saddleNum*2);
	integrationBackwardCurves.setParameters(0.5, 2000);
	backward = integrationBackwardCurves.integrate();
}

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
/* Unfortunately, the algorithm used to find zero point is unfatisfactory.
 * This is because the zero flow is a necessary but not sufficient
 * condition that the point is a critical point.
 */
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
bool visFlowTopology::isZeroPoint(size_t id, double *pos)
{
	vector<double> dataV;
	dataField->getDataValueOfPoint(id, dataV);
	int m = 0;
	for(; m < 2; m++)
	{
		if( abs(dataV[m]) > 1e-6)
		{
			break;
		}
	}
	// If m >= 2 holds, it can be known that the vector dataV is zero.
	if(m >= 2)
	{
		vector<double> coords;
		dataField->getPointCoordinates(id, coords);
		for(int i = 0; i < 3; i++)
		{
			pos[i] = coords[i];
		}
		return true;
	}
		
	/* For each point thita on the circle
	 *   1. calculate the local coordinate of thita.
	 *   2. translate it to the global coordinate by adding the coordinates
	 *      of thita to that of the point id.
	 *   3. interpolate to get the vector pntOnCircle of point thita.
	 *   4. calculate the dot product of the two vectors pntOnCircle and
	 *   	(id -> thita).
	 *   5. add the result to variable flow.
	 */
	double angle = 2 * 3.14156/smpRate;

	vector<double> pointCoordsOfId;
	dataField->getPointCoordinates(id, pointCoordsOfId);

	double pntOnCircle[][3] = 
	{
		{0, 0, 0},
		{0, 0, 0},
		{0, 0, 0},
		{0, 0, 0}
	};

	int i = 0;
	vector<double> dataVOfPnt[4];
	for(double thita = 0; thita < 2 * 3.14156; thita += angle)
	{
		double consine = cos(thita);
		double sine = sin(thita);
		pntOnCircle[i][0] = pointCoordsOfId[0] + radius * consine;
		pntOnCircle[i][1] = pointCoordsOfId[1] + radius * sine;

		// For the data set 19, the nearest number should be 5.
		//if(!dataField->interpolatedValue(pntOnCircle[i], 3, &dataVOfPnt[i], 4, 1.5))
		wunai = 4;
		if(!dataField->interpolatedValue(pntOnCircle[i], 3, &dataVOfPnt[i], wunai, 4))
		{		return false; 	}
		i++;
	}

	/*           1
	 *           |
	 *      2----c----0
	 *           |
	 *           3  
	 *
	 */
	if( ( dataVOfPnt[0][0] * dataVOfPnt[2][0] + // 0, 2 : x
	      dataVOfPnt[0][1] * dataVOfPnt[2][1] < 0 // 0, 2 : y
	     )
	  || 
	    ( dataVOfPnt[1][0] * dataVOfPnt[3][0] + // 1, 3 : x
	      dataVOfPnt[1][1] * dataVOfPnt[3][1] < 0 // 1, 3 : x
	    )
	   )
	{
		for(int i = 0; i < 3; i++)
		{	pos[i] = pointCoordsOfId[i];	}

		return true;
	}

	return false;
}

bool visFlowTopology::computeJacobian(double *pos)
{
	/* 1. Get coordinates of point whose id is "id";
	 * 2. Calculate the differences of u and v along x and y respectively. That is,
	 *    ux, uy, vx, vy.
	 * The matrix is :
	 *  /           \
	 *  | ux    uy  |
	 *  | vx    vy  |
	 *  \           /
	 *  or, [ ux, vx, uy, vy].
	*/

	// 1.
/*	double step = delt;
	std::vector<double> coordinates;
	dataField->getPointCoordinates(id, coordinates);

	double pos[3];
	for(unsigned i = 0; i < coordinates.size(); i++)
	{
		pos[i] = coordinates[i];
	}
*/
	// 2.
	vector<double> uvXPlus;
	vector<double> uvXSub;
	vector<double> uvYPlus;
	vector<double> uvYSub;	

	pos[0] += delt;
	if(!dataField->interpolatedValue(pos , 3, &uvXPlus, neighbourNum, disError))
	{
		return 0;
	}
	pos[0] -= delt;

	pos[0] -= delt;
	if(!dataField->interpolatedValue(pos , 3, &uvXSub, neighbourNum, disError))
	{
		return 0;
	}
	pos[0] += delt;

	pos[1] += delt;	
	if(!dataField->interpolatedValue(pos , 3, &uvYPlus, neighbourNum, disError))
	{
		return 0;
	}
	pos[1] -= delt;

	pos[1] -= delt;
	if(!dataField->interpolatedValue(pos , 3, &uvYSub, neighbourNum, disError))
	{
		return 0;
	}
	pos[1] += delt;

	jacobian[0] = (uvXPlus[0] - uvXSub[0]) / (2 * delt);
	jacobian[1] = (uvXPlus[1] - uvXSub[1]) / (2 * delt);
	
	jacobian[2] = (uvYPlus[0] - uvYSub[0]) / (2 * delt);	
	jacobian[3] = (uvYPlus[1] - uvYSub[1]) / (2 * delt);

	return 1;
}

bool visFlowTopology::matrixIsZero(const double jac[dim * dim])
{
	for(int i = 0; i < dim * dim; i++)
	{
		if(abs(jac[i]) > zero)
		{
			return 0;
		}
	}
	return 1;
}

//bool visFlowTopology::calculator(double jac[], complex_ eigval[], complex_ eigvec[])
//{
//	// If dim is not 2, the program considers an error comes, since I'am now
//	// processing two dimensional cases.
//	if(2 != dim)
//	{
//		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
//		cout << "The dimension should be 2, but not others ..." << endl;
//		exit(-1);
//	}
//	integer n = dim;
//	complex_ a[dim * dim];
//	for(int i = 0; i < dim * dim; i++)
//	{
//		a[i].r = jac[i];
//		a[i].i = 0;
//	}
//
//	integer lda = dim;
//	complex_ vl[dim * dim];
//	integer ldvl = dim, ldvr = dim;
//	complex_ work[2 * dim];
//	real_ rwork[2 * dim];
//	integer info;
//	integer lwork = 2 * dim;
//
//	cgeev_((char*)"N", (char*)"V", &n, a, &lda, eigval, vl, &ldvl, eigvec, &ldvr, work, &lwork, rwork, &info);
//	if(0 != info)
//	{
//		cout << "AT LINE " << __LINE__ << " OF FILE " << __FILE__ << endl;
//		cout << "EIGENVALUE AND EIGENVECTOR FAILS!" << endl;
//		return 0;
//	}
//// The following code is to display the calculated eigenvectors and eigenvalues.
///*
//	cout << "The eigenvalue: " << endl;
//	for(int i = 0; i < dim; i++)
//	{
//		cout << eigval[i].r << " + " << eigval[i].i << "i" << endl;
//	}
//
//	cout << "The eigenvector: " << endl;
//	for(int i = 0; i < dim * dim;)
//	{
//		cout << "[ ";
//		for(int j = i; j < i + dim; j++)
//		{	cout << eigvec[j].r << " + " << eigvec[j].i << "i";
//			if(j - i < dim - 1)
//			{
//				cout << ", ";
//			}
//		}
//		cout << " ]" << endl;
//		i += dim;
//	}
//*/
//	return 1;
//}

bool visFlowTopology::classify(const double *pos)
{
 //Type{ATTRACTINFOCUS, REPELLINGFOCUS, ATTRACTINGNODE, REPELLINGNODE, CENTER};

	const bool isSaddle = 1;
	// If the imaginary parts of both of the eigenvalues are non-zero.
	if(zero < abs(eigenvalue[0].imag()) && zero < abs(eigenvalue[1].imag()))
	{
		// repelling focus
		if(eigenvalue[0].real() > 0 && eigenvalue[1].real() > 0)
		{
			push(criticalPoints, pos, REPELLINGFOCUS);
			return !isSaddle;
		}
		// attracting focus
		if(eigenvalue[0].real() < 0 && eigenvalue[1].real() < 0)
		{
			push(criticalPoints, pos, ATTRACTINFOCUS);
			return !isSaddle;
		}
		// center
		if(abs(eigenvalue[0].real()) < zero && abs(eigenvalue[1].real()) < zero)
		{
			push(criticalPoints, pos, CENTER);
			return !isSaddle;
		}
	}
	else if(abs(eigenvalue[0].imag()) < zero && abs(eigenvalue[1].imag()) < zero)
	{
		// saddle!
		if(eigenvalue[0].real() * eigenvalue[1].real() < 0)
		{
			saddlePoint tmp;
			push(tmp, pos);
			
			saddles.push_back(tmp);
			return isSaddle;
		}
		// repelling node
		if(eigenvalue[0].real() > 0 && eigenvalue[1].real() > 0)
		{
			push(criticalPoints, pos, REPELLINGNODE);
			return !isSaddle;
		}
		// attracting node
		if(eigenvalue[0].real() < 0 && eigenvalue[1].real() < 0)
		{
			push(criticalPoints, pos, ATTRACTINGNODE);
			return !isSaddle;
		}
	}
	// None of the above is satisfies. Note that this is possible.
	return 0;
}

void visFlowTopology::normalizeVec(double vec[dim * dim])
{
	double sum;
	sum = pow(vec[0], 2) + pow(vec[1], 2);
	// Note that the "sum" can be zero. If so, it will be unnecessary
	// to normalize the corresponding vector.
	if(1e-15 < sum)
	{
		vec[0] /= sqrt(sum);
		vec[1] /= sqrt(sum);
	}

	sum = pow(vec[2], 2) + pow(vec[3], 2);
	if(1e-15 < sum)
	{
		vec[2] /= sqrt(sum);
		vec[3] /= sqrt(sum);
	}
}

}// visualization

